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Introduction 

The  ambiguity  function  or  surfaces  provides  the  designer  of  modern  radar  and  sonar 
technology  a  measure  of  the  effectiveness  of  a  given  matched-filter  signal  in  determining  the 
simultaneous  evaluation  of  range  and  velocity  of  a  moving  object  to  the  particular  conditions 
of  target  environment.  In  this  work,  we  will  be  concerned  with  efficient  computational 
methods  for  evaluating  this  function  by  computers. 

The  ambiguity  function  is  defined,  for  our  purposes,  as  the  absolute  value  of  the  function 

G(j,  n)  -  5)  xky^+j  exp  (-2tri  -2^) 

k-o  K 

where  J,  J£,  N  are  positive  integers,  0<k<K,  0<j<J,  -N<n<N  and  xk  yk+j  are  data.  We 
will  denote  1  by  i. 

In  section  2.,  we  will  view  the  computation  of  G(j,n)  as  a  filter  or  convolution 
nk  * 

(of  xk  exp  (~2*ri— )  by  yk),  while  in  section  3.,  G(j,n)  is  interpreted  as  a  Fourier  transform 
(of  xk  yk+1).  Several  departures,  from  usual  methods,  making  use  of  special  circumstances, 
will  be  described,  significantly  increasing  the  computational  efficiency. 

In  section  3.,  as  well,  the  technique  of  passing  a  long  sequence  through  an  F.I.R.  filter 
with  decimation  and  computing  the  discrete  Fourier  transform  (D.F.T.)  of  the  shorter  decimat¬ 
ed  sequence  will  be  applied  to  the  problem.  We  will  follow  the  approach  of  [  1  ]. 

In  section  4.,  a  method  based  on  polynomial  approximation  theory  will  be  introduced  and 
compared  to  the  decimating  filter  technique.  Indeed,  after  some  initial  differences,  the  final 
stages  of  computing  by  the  approximation  technique  coincides  with  those  of  the  decimating 


filter  method. 
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2.  Filter  Method 

Let  A(j,n)  ■  |  G(j,n)  | .  The  first  method  for  computing  A(j,n)  is  based  upon  writing 
G(j,n)  as  a  filter  and  applying  fast  Fourier  transform  (FFT)  techniques  to  the  resulting 
computation.  Let 


zk(n)  —  xk  exp  (  -  2wi— ). 

1C 

For  each  n,  -N<n<N,  we  can  write 

K-l 

(2.1)  G(j,n)  «  2  zk(n)  yk+j  ,  0<j<J 

k-0 


which  expresses  G(j,n)  as  the  first  J+l  outputs  of  a  K-tap  filter.  (That  is,  we  view  the 
sequence  {y^}  as  the  data  and  the  sequence  {zk(n>}  as  the  tap  values.) 

Equivalently,  if  we  set 


,0£k<K 
,  K<k<K  +  J 


(padding) 


^  • 

y ;  *  Yk+j-/  ,0<^<K-*-J  (bit -reversal) 


then  for  each  n,  —  N<n<N,  we  can  write  G(j,n)  as  the  cyclic  convolution. 
(2  2)  G(j,n)  -  (z(n)  •  y)K+J_j  ,  0<j<J. 


Computing  the  cyclic  convolution  by  standard  DFT  techniques  yields  Flow  chart  A  for 
the  computing  of  G(j,n),  0<j£ J.  In  general,  the  DFT  of  a  sequence  x  will  be  denoted  by 
X. 


1 
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Denote  by  FFT(K)  the  number  of  complex  multiplications  used  by  the  FFT  methods  to 

~  • 

compute  the  DFT(K).  Note  y  is  independent  of  n  and  Y  need  be  computed  only  once  in 
Flow  chart  A.  It  follows  that  the  number  of  complex  multiplications  used  by  Flow  chart  A  to 
compute  G(j,n),  0<j<J,  -N<n<N,  is  given  by 

(2.3)  (2N  +  1)  (K  +  2  FFT(K  +  J)  +  (K  +  J))  +  FFT(K  +  J). 

The  number  given  by  (2.3)  will  be  denoted  by  w,  —  w,(N,  K,  J).  Table  I  describes  the 
numerology  for  three  choices  of  N,  K  and  J.  In  the  examples,  we  will  take 
FFT(K)  -  Ik  log2  K. 

Table  I 


K 

J 

N 

n, 

213 

2*3 

50 

3145  .  213 

2U-700 

700 

50 

1623  .  214— 70,700  ~  1619  .  2U 

2U-70 

70 

50 

1623  •  214— 7,070  ~  1623  .  214 

In  certain  cases,  we  can  improve  on  the  computational  efficiency  of  the  previous  method. 
Suppose,  for  instance,  we  have  K  -  J.  The  computation  of  Z(n),  -N<n<N,  can  be  made  as 
follows. 


Pad  the  sequence  x  by  setting 


(xk  ,  0<k<K 
0  ,K£k<2K 
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and  denote  the  DFT(2K)  of  x  by  X.  Then, 

X„+2„  -  x,  exp(-2irii^k) 

k-0 

-  j  zk<“>  exp(-2iri^-) 
k-0 

2K-1  _ 

-  2  zk<n)  exp(-2wi^). 

k-0  ZK 

It  follows  that 

<v  M 

(2.4)  (Z(n))u  -  Xu+2n  ,  0<u<2K,  — N<n<N. 

Once  X  has  been  computed,  without  any  additional  complex  multiplications,  we  can 
compute  Z(n),  -N<n<N.  Specifically,  Z(n)  is  given  by  cyclically  shifted,  mod  2K,  the 
sequence  X  2n  places. 

Flow  chart  B  descibes  how  the  observation  applies  to  the  computation  of  G(j,n).  Also, 
the  number  of  complex  multiplications  used  on  Flow  chart  B  is  given  by 

(2.5)  (2N  +  1)(2K  +  FFT(2K))  +  2  FFT(2K). 
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Example.  For  K  ■  J  m  213  and  N  ■  50,  substituting  into  (2.5)  yields  1644.2 13  which  is  a 
50%  savings  as  compared  to  ir,  as  given  in  Table  I. 

! 

Essentially,  the  same  argument  holds  if  J  is  any  integer  multiple  of  K.  The  case  when  K 
ia  a  multiple  of  J  will  now  be  considered.  Suppose,  for  instance,  K  -  2J.  Pad  x  by  setting 


(tk  ,  0<k<K 
(0  ,K<k<lK 


and  let  X  be  the  DFT(—  K)  of  x,  as  above, 

2 

(2.6)  (Z(2n))u  .  Xu+3n  ,  0<u<-|  K,  -N<2n<N. 

it  follows  that  Z(2n)  can  be  found  by  cyclically  shifting,  mod  (—  K),  the  sequence  X  3n 

2 

places.  Analogously,  if  we  set 


exp(-2wi 

0 


,  0<k<K 

,  K<k<— K 
2 


then 


(2.7)  (Z(2n  +  l))u  -  X(l)u+3n 

Thus,  the  number  of  complex  multiplications  used  is 


(2.8)  (2N  +  1)(FFT(—  K)  +  —  K)  +  3  FFT(—  K)  +  K. 

2  2  2 


Example  2.  Suppose  K  ■  213,  J  ■  2 12  and  N  ■  50.  Flow  chart  A  uses  approximately 
2>z  .  6000  complex  multiplications.  By  (2.7),  since  K  ■  2J,  we  can  compute  the  case  in 
ar  ximately  2 12 . 4650  complex  multiplications. 
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3.  Transform  Method 

For  each  j,  0<j<J,  consider  the  sequences,  defined  mod  K, 


yk(j)  -  yk+j 


uk(j)  -  *kyk(j)- 


,  0<k<K, 


We  can  write 

(3.1)  G(j,  n)  -  U„(j)  ,  -N<n<N, 

where  U(j)  is  the  DFT(K)  of  u(j).  Flow  chart  C  outlines  the  computation  of  G(j,n)  by  this 
method.  The  number  of  complex  multiplications  used  by  the  method  is 

(3.2)  (J  +  1)(FFT(K)  +  K). 

We  will  denote  the  number  given  by  (3.2)  by  rr2  m  W;(K,  J).  Table  II  reconsiders  those 
examples  of  Table  I  by  this  new  method  and  compares  these  values  corresponding  values  ir, 
given  in  Table  I. 


Flow  Chart  C 


x  y 


y(j) 


u(j) 


DFT(K) 


U(j) 

Table  II 


K 

J 

N 

w2 

ffi 

2'3 

213 

50 

212  .  15  .  8193 

214— 700 

700 

50 

<  5586  .  214 

2l4-70 

70 

50 

<  567  .  214 

~  3 

Note  that  in  the  last  case  the  transform  method  is  more  computationally  efficient  than  the 
filter  method. 

In  the  most  important  applications,  N<<K  and  we  will  assume  this  in  all  that  follows. 
As  evidenced  by  (3.2),  the  method  of  Flow  chart  C  does  not  take  into  account  the  relative 
size  of  N  as  compared  to  K.  Generally  speaking,  the  problem  is  to  compute  the  DFT(K)  of 
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the  long  sequence  u(j)  for  only  a  small  interval  of  frequencies  — N<n<N.  We  will  follow  the 
approach  of  [1].  As  we  will  see,  the  long  sequence  u(j)  will  be  passed  into  an  FIR  filter  with 
decimation  and  the  DFT  of  the  shorter  decimated  sequence  will  be  computed.  The  amount  of 
complex  multiplications  used  will  depend  upon  the  number  of  taps  of  the  filter  and  the 
decimation  rate.  However,  the  DFT  of  the  decimated  sequence  yields  only  an  approximation 
to  G(j,n),  due  to  aliasing  errors. 

Consider  to  begin  with,  a  K-tap  filter  given  by 

h  -  (h(0),  h(l) . h(K-l)) 

satisfying  the  following  properties 

(I)  h(k)  -  h(-k). 

(II)  h(k)  vanishes  outside  an  interval,  -M<k<M,  for  some  integer  M. 

(III)  Let  H  be  the  DFT(K)  of  h.  Then  H(0)  -  1  and  H(n)  never  vanishes,  -N<n<N. 

Passing  the  sequences  u(j)  through  the  filter  b  we  determine  the  sequence  v(j)  given  by 

K-l 

v,(j)  -  2  Wi>h<k>  ,  0<f <K, 

k-0 

or  equivalently  v(j)  *  u(j)*h.  Let  U(j),  V(j)  denote  the  DFT(K)  of  u(j),  v(j),  respectively. 
Then, 

V(j)  -  U(j)  x  H. 


Suppose  K  ■  R  x  L,  R,  L  integers  and  assume  N<R.  Consider  the  (L:l)  decimation  of 


br(j)  -  vrL(j) 


0<r<R 


and  denote  the  DFT(R)  of  b(j)  by  B(j).  A  standard  computation,  see  [1],  shows  that 

Br(j)  -  L_,(Vr(j)  +  Er(j))  ,  0<r<R, 


where 

L—  1 

Er(j)  -  2  Ur+*R<i>  H(r  +  *R)  ,  0<r<R. 

4-1 

Condition  (III)  implies  we  can  write 

Un(j)  -  H(n)-1  L  Bn(j)— H(n)-1  E„(j)  ,  -N<n<N. 

We  plan  to  approximate  G(j,n)  ■  Un(j)  by  H(n)-IL  Bn(j).  The  error 

i  L“‘  i 

H(n)  ' E„(j)  -  2  Un+4R(i>  H(">  H<n  +  *R> 

4-1 

is  due  to  the  aliasing  of  Un+4R(j)  into  the  frequencies  — N<n<N  attenuated  by  the  factors 
H(n)  1  H(n  +  4R),  — N<n<N,  0<4<L.  An  analysis  of  this  aliasing  error  was  carried  out  in 
[1].  For  our  purposes,  it  suffices  to  remark  that  the  filter  h  should  be  chosen  to  minimize 
the  factors  H(n)  1  H(n  +  4R),  for  -N<n$N,  0<4<L,  in  the  sense  of  the  sup  norm. 
Condition  (I)  permits  the  application  of  the  Remez  algorithm  as  contained  in  [1]  to  achieve 
the  minimization  of  these  factors  subject  to  conditions  (II)  and  (III).  At  the  same  time,  the 
Remez  algorithm  produces  the  uniform  bound  of  these  factors. 

For  each  M  and  L,  the  Remez  algorithm  determines  a  filter  h  satisfying  (I),  (II)  and 
(III)  and  minimizing  the  factors  H(n)~’  H(n  +  4R)  over  the  aliasing  frequencies.  For  fixed  L. 
the  error  term  increases  as  the  number  of  taps,  2M  +  1,  of  h  decreases.  Also,  increasing 
the  decimation  rate  L,  with  fixed  M,  increases  the  error  term.  Thus,  M  large  and  L  small 


produces  small  error  terms.  However,  as  we  will  soon  see,  the  computation  of  H(n)-IL  Bn(j) 
is  most  computationally  efficient  when  M  is  small  and  L  is  large. 

Variations  of  these  ideas  are  considered  in  [1]  where  weighting  factors  over  the  decreas¬ 
ing  frequencies  are  introduced  into  the  Remez  subroutine  as  a  way  of  accounting  for  the  a 
priori  knowledge  that  certain  aliasing  frequencies  require  less  attenuation  than  others. 

Flow  chart  D  describes  the  computation  of  H(n)  *L  Bn(j),  in  a  straightforward  manner. 
A  more  computationally  efficient  method  will  eventually  be  given. 
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Flow  Chart  D  ForO<j<J, 

x  y*(j) 


H(n)-1  L  B„(j). 


The  number  of  complex  multiplications  used  by  Flow  chart  D  to  compute 
H(n)-1  L  Bn(j),  -NSngN,  0£j $2J  is  given  by 


(3.3) 


(J  +  1)(K  +  2M  +  1  +  FFT(R)  +  2N  +  1) 
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The  major  portion  of  the  computation  occurs  in  passing  from  x  and  y(j)  to  B(j).  We  will 
study  this  process  more  closely.  First,  we  can  write 

K-l 

br(j)  -  2  urL-k(i>  h<k> 

k-0 

K-l 

"  2  XrL-k  yrL-k(j) 

tm  0 

M 

2  xrL-rnh<m>  y^L-mO) 

m— — M 

Set  fm(r)  -  xrL+M_mh(  -  M  +  m),  0$m<2M  and  gm(r)  «  y*L_M+m,  0<m<2M  +  J. 
Then, 

2M 

(3-4)  br(j)  -  £  fB(r)  gm+j(r). 

m>0 

The  process  described  in  Flow  chart  D  assumes  j,  0<j$J  is  fixed.  We  will  now  adopt  the 
point  of  view  that  r,  0£r<R,  is  fixed  and  use  (3.4)  to  compute  br(  ),  the  sequence  in  j, 
0<j<J,  corresponding  to  the  fixed  r.  Thus,  (3.4),  expresses  br(j)  as  the  first  (J+l)-outputs 
of  a  (2M+l)-tap  filter  and  the  methods  of  section  2.,  especially  Flow  chart  A,  can  be 
employed.  We  are  also  assuming  M  is  sufficiently  large  to  make  the  approach  computationally 
efficient. 

Flow  chart  E,  below,  describes  how  for  fixed  r,  0<r<R,  the  method  of  Flow  chart  A 
proceeds  from  original  data  to  the  sequence  br(  ).  If  Flow  chart  E  is  carried  out  for  each  r, 
0<r<R,  then  the  end  of  Flow  chart  D  can  be  applied  to  b(j)  to  compute  H(n)-IL  Bn(j). 


r 


l 


DFT(2M+1+J) 


DFT(2M+1+J) 
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Flow  chart  E  uses 

(3.5)  R(2M  +1  +  3  FFT(2M  +  1  +  J)  +  2M  +  1  +  J) 

complex  multiplications  to  compute  bf(j),  0$j<J,  0<r<R,  and  hence 
H(n)~*L  Bn(j),  -N$n<N,  0<j$J,  can  be  computed  by 

(3.5)  R(2M  +1  +  3  FFT(2M  +  1  +  J)  +  2M  +  1  +  J)  +  (J  +  1)(FFT(R)  +  2N  +  1) 

complex  multiplications. 

In  Table  III  below,  we  will  denote  the  number  given  by  (3.6)  by  v,  -  w3(N,  K,  J,  M,  L). 

Table  3  N  -  50 

K  J  M  L  ir3  w2 

2*3  2*3  29-l  27  >  212  •  508  212  .  15  .  8193 

For  a  more  complete  analysis  of  the  relationship  between  given  values  of  K.  J,  M,  L,  N 
and  the  corresponding  values  of  ir3  see  Table  I,  pages  216  of  [1]. 

4.  Approximation  Theory 

In  this  section,  we  will  obtain  a  method  analogous  to  that  discussed  in  section  3,  by  using 
ideas  from  approximation  theory.  As  in  section  3,  the  DFT(R)  will  be  applied  to  the  (L:  1 ) 
decimation  of  a  sequence  obtained  by  passing  K-point  data  through  a  filter.  The  filter  will  be 
determined  by  the  degree  of  the  approximation  taken,  in  a  sense  explained  below.  Once  the 
computation  to  be  made  is  written  in  this  form,  the  methods  of  section  3  can  be  employed. 


(4.5) 


(J  +  l)(K  +  FFT(R)). 


Clearly,  (4.4),  expresses  cr(j)  as  the  (L:l)  decimation  of  the  filtered  sequence 


wkW  *  2  “k+fW  ,0<k<K. 

tm  o 


For  n,  —  N<n<N,  consider  that  part  of  the  unit  circle  given  by 

f(t)  m  exp(-2iri— ,  0<t<L.  Denote  by  (t)  the  polynomial  of  degree  one  in  t  uniquely 
1C 

determined  by  the  conditions 


q,(0)  -  f(0) 
qi(L)  -  f(L). 


Thus,  q,(t)  is  geometrically  the  line  joining  the  point  1  to  the  point  exp(-2iri  — ).  We  can 

K 


~2*i-2- 


q,(t)  -  (1-£>  +  £e  * 


Choose  a  real  number  d,(n)  such  that 


P|(t)  -  dj(n)  q,(t) 


is  the  best-approximation  to  f(t),  in  the  sup-norm  sense,  over  the  interval  0<t<L,  within  the 


space  of  ail  real  multiples  of  gj(t). 


Replacing  exp(-2iri-^-)  by  p , (/)  in  (4.2)  we  obtain 


d,(n)  Cn(j) 


where 


-  2  (urL+f(i>(,-f)  +  U(r-DL+^>{) 
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and  C(j)  is  the  DFT(R)  of  c(j). 

To  see  the  analogy  of  this  approach  to  the  method  of  section  3,  we  proceed  as  follows. 
Let  h(k),  0<k<K,  be  defined  by 

,0£k<L 

±  ,  L<k<2L 

L 

and  zero  otherwise.  Consider  the  2L-tap  filter  determined  by  h, 

2L-1 

(4  8)  *  wt(j)  *  2  M*>  u,_L+/(i) 

r-0 

By  (4.7)Twe  can  write  c(j)  as  the  (L:l)  decimating  filter 
(4.9)  cr(j)  -  wrL(j)  ,  0<r<R. 

Again,  we  are  led  to  the  computing  of  the  DFT(R)  of  an  (L:l)  decimating  filter.  Here,  the 
filter  is  determined  by  the  approximation  taken.  Since  the  interpolating  polynomial  has  degree 
one,  we  call  the  method,  the  first  degree  approximation  method. 

The  techniques  of  section  3,  can  be  applied  to  the_computation  of  cr(j).  Setting 

X/(r)  -  x(r_i)L+/ 

x',(r)  -  x,(r)  h(0 

y'j +/(f)  -  y](r)  y'/+j(r)’ 

which  for  each  r,  0£r<R,  expresses  cr(  )  as  the  first  (J4-l)-outputs  of  a  2L-tap  filter.  The 
number  of  complex  multiplications  used  by  this  method  to  compute 
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d,(n)  C„(j),  -N<n<N,  0<j<J,  is 

(4.11)  2K  +  R(3FFT(2L  +  J)  +  2L  +  J)  +  (J  +  1)(FFT(R)  +  2N  +  1). 

Choose  M>0  and  consider  that  part  of  the  unit  circle  given  by 

f(t)  m  exp(-2vi  — ),  0<t<ML.  Denote  by  qM(t)  the  uniquely  determined  M-th  degree 
1C 

polynomial  in  t  satisfying 

qM(0  -  f(t)  ,  t  -  0,  L,...,  ML. 

Using  Lagrange  interpolation  we  can  write 

W1)  “  2  hm(t)  exp(-2iri  Z&-) 

m-0  K 

where  h0(t) . hM(t)  are  polynomials  in  t.  Choose  a  real  number  dM(n)  such  that 

pM(t)  ■  dM(n)  qM(t) 

is  the  best-approximation  to  f(t),  in  the  sup-norm  sense,  over  the  interv-1  0<t<ML,  within 
the  set  of  real  multiples  of  qM(t).  We  call  pM(t)  an  M-th  degree  approximation  since  the 
interpolating  polynomial  has  degree  M. 

Substituting  pM(0  for  exp(-2tri  in  (4.2),  yields 

(4.12)  dM(n)  CB(j) 

where 

^  13)  cr<j>  -  2  2  hm<'> 

fmO  m-0 


and  C(j)  is  the  DFT(R)  of  c(j). 


The  analogy  with  section  3  can  be  seen  by  defining  h(k),  0<k<K,  by 


Arguing  as  in  the  previous  case,  we  can  express  the  sequence  (cr(0) . cr(J)>  as  the  first 

(J+  l)-outputs  of  an  (M+l)  •  L-tap  filter.  The  number  of  complex  multiplications  used  to 
compute  dM(n)  Cn(j),  -N<n^N,  0<j£j,  by  this  method  is 

(4.15)  (M  +  1)K  +  R(J  +  (M  +  1)L  +  3FFT(J  +  (M  +  1)L))  +  (J  +  1)(FFT(R)  +  2N  +  1). 
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